# Importamos librerias necesarias para nuestro programa
import pandas as pd
import numpy as np
import plotly.express as px
import plotly.graph_objects as go  # Ensure plotly.graph_objects is imported
from IPython.display import display, Markdown
from scipy.signal import find_peaks
from sklearn.metrics import r2_score
from sklearn.linear_model import LinearRegression
from scipy.optimize import curve_fit
import plotly.io as pio
pio.renderers.default = "png"

file = open('../SrTiO3.uxd', mode='r')

content = file.read()

partes_importantes = content.split(';')

tabla_contenido = partes_importantes[7]

titles = tabla_contenido[1:15]

tabla_contenido = tabla_contenido.replace(titles, '2THETA, PSD\n')
tabla_contenido = tabla_contenido.replace('       ', ', ')
tabla_contenido = tabla_contenido.replace('      ', ', ')

file.close()

output_file = open('data.csv', 'w')

output_file.write(tabla_contenido)

output_file.close()

df = pd.read_csv('data.csv')

y_position_value_global = 100
beta_constant_value = 0


fig = px.line(df, x=' 2THETA', y=' PSD', labels={'Name': '2theta', 'Value': 'values'}, title='SrTiO3 Difractograma')
fig.update_traces(line=dict(color='blue')) 

# Add a line parallel to x-axis at y_position_value
fig.add_shape(
    type='line',
    x0=df[' 2THETA'].min(),
    y0=y_position_value_global,
    x1=df[' 2THETA'].max(),
    y1=y_position_value_global,
    line=dict(color='red', width=2, dash='solid')
)

#-------------------------------------------------
#Se calcula la altura máxima en de la reflexión más grande:
altura_reflexion_max = max(df[' PSD'].values) - y_position_value_global
# indice reflexión más alta
indice_refle_mas_alta = df[df[' PSD'] == max(df[' PSD'].values)].index[0]
#print(indice_refle_mas_alta)
# Obtenemos la mitad de la distancia con de la reflexión mas grande usando el punto de referencia
given_PSD_value = altura_reflexion_max/2

# Se calcula la diferencia absoluta entre los dos valores dados y todos los datos de la columna PSD.
df['Absolute_Difference'] = abs(df[' PSD'] - given_PSD_value)


# Encuentra la fila con la menor diferencia encontrada
closest_index = df['Absolute_Difference'].idxmin()

closest_2THETA_value = df.loc[closest_index, ' 2THETA']
closest_PSD_value = df.loc[closest_index, ' PSD']

# Encontramos el siguiente valor más cercano
next_upper_values = df[df[' PSD'] > given_PSD_value]
if not next_upper_values.empty:
    next_upper_index = next_upper_values[' PSD'].idxmin()
    next_upper_2THETA = df.loc[next_upper_index, ' 2THETA']
    next_upper_PSD = df.loc[next_upper_index, ' PSD']

#-------------------------------------------------

promedio_de_dos_puntos = ((df[' 2THETA'][next_upper_index] + df[' 2THETA'][next_upper_index + 1])/2) + 0.003

# Agregamos dos puntos con las coordenadas de closes index y el promedio de 2theta en el next_upper_index y next_upper_index + 1
fig.add_trace(go.Scatter(x=[df[' 2THETA'][closest_index]], y=[df[' PSD'][closest_index]], mode='markers+text', text='Punto A', textposition='bottom center', marker=dict(color='blue')))
fig.add_trace(go.Scatter(x=[promedio_de_dos_puntos], y=[df[' PSD'][closest_index]], mode='markers+text', text='Punto B', textposition='bottom center', marker=dict(color='red')))


# Creamos una linea entre esos cos puntos
fig.add_trace(go.Scatter(x=[df[' 2THETA'][closest_index], promedio_de_dos_puntos], y=[df[' PSD'][closest_index], df[' PSD'][closest_index]], mode='lines', line=dict(color='green', dash='dash')))


# Se calcula la distancia con la formula euclidiana
def calculate_fwhm(index):
    half_max = df[' PSD'][index] / 2  # Half of the peak's maximum value
    left_bound = np.where(df[' PSD'][:index] < half_max)[0][-1]  # Left boundary
    right_bound = np.where(df[' PSD'][index:] < half_max)[0][0] + index  # Right boundary
    fwhm = df[' 2THETA'][right_bound] - df[' 2THETA'][left_bound]  # FWHM calculation
    return fwhm

beta_constant_value = calculate_fwhm(indice_refle_mas_alta)  #Se guardan los datos en la variable global
distance_two_points = beta_constant_value
# Se agrega una anotación con esta distancia
fig.add_annotation(
    x=(next_upper_2THETA + closest_2THETA_value)/2,
    y=df[' PSD'][closest_index],
    text=f'Distancia: {distance_two_points:.5f}', # Ponemos la distancia con 3 puntos decimales
    showarrow=True,
    arrowhead=1,
)

fig.add_annotation(
    x=(next_upper_2THETA + closest_2THETA_value)/2,
    y=df[' PSD'][closest_index + 3],
    text=f'Distancia entre linea y pico: {altura_reflexion_max} y distancia mitad {altura_reflexion_max/2}', # Ponemos la distancia con 3 puntos decimales
    showarrow=True,
    arrowhead=1,
)

fig.show(renderer='notebook')

def calc_tamanio_cristalito_scherrer (beta, theta):
    K = 0.89 # Para cúbicas según la referencia
    LAMBDA = 1.5406 # Longitud de onda Cobre K alfa
    
    #print(theta)
    theta_rads = np.deg2rad(theta)
    angulo = np.cos(theta_rads)
    #print(angulo)
    cristalito_size = (K * LAMBDA)/(beta * angulo)
    #print(f'Tamaño de cristalito calculado: {cristalito_size.values[0]} Å')
    return cristalito_size.values[0] * 0.1 # Se multiplica por 0.1 para convertir A a nm

theta_scherrer = df[df[' PSD'] == max(df[' PSD'].values)][' 2THETA']
#print(theta_scherrer/2)
display(Markdown(f'''
&nbsp;

<div align="center"> Tamaño de cristalito calculado: <b> {calc_tamanio_cristalito_scherrer(np.deg2rad(distance_two_points),(theta_scherrer/2))} nm </b></div>

&nbsp;

'''))

# Encontramos los picos más altos en la gráfica
peaks, _ = find_peaks(df[' PSD'], prominence=50 , distance= 90)  # Adjust prominence as needed
peaks = np.delete(peaks, 0)
peaks = np.delete(peaks, 1)
peaks = np.delete(peaks, 7)

y_position_value_global = 100
beta_constant_value = 0

lista_picos_index = []

for i in peaks:
    lista_picos_index.append(i)

fig = px.line(df, x=' 2THETA', y=' PSD', labels={'Name': '2theta', 'Value': 'values'}, title='SrTiO3 difractograma')
fig.add_scatter(x=df[' 2THETA'][peaks], y=df[' PSD'][peaks], mode='markers', marker=dict(color='red', size=8), name='Peaks')

fig.update_traces(line=dict(color='blue')) 



# Agregamos una linea horizontal que sirve como punto de referencia
fig.add_shape(
    type='line',
    x0=df[' 2THETA'].min(),
    y0=y_position_value_global,
    x1=df[' 2THETA'].max(),
    y1=y_position_value_global,
    line=dict(color='red', width=2, dash='solid')
)

#-------------------------------------------------
#Se calcula la altura máxima en de la reflexión más grande:
#altura_reflexion_max = max(df[' PSD'].values) - y_position_value_global


def calculate_fwhm(index):
    half_max = df[' PSD'][index] / 2  # Half of the peak's maximum value
    left_bound = np.where(df[' PSD'][:index] < half_max)[0][-1]  # Left boundary
    right_bound = np.where(df[' PSD'][index:] < half_max)[0][0] + index  # Right boundary
    fwhm = df[' 2THETA'][right_bound] - df[' 2THETA'][left_bound]  # FWHM calculation
    return fwhm

def largo_pico(distancia, df_f, peaks):
    
    fig.add_annotation(
        x=df_f[' 2THETA'][peaks],
        y=df_f[' PSD'][peaks - 2],
        text=f'FWHM {distancia:.3f}', # Ponemos la distancia con 3 puntos decimales
        showarrow=True,
        arrowhead=1,
    )

fwhm_values = []
angulos_peaks = []
for peak_index in peaks:
    fwhm = calculate_fwhm(peak_index)
    angulos_peaks.append(df[' 2THETA'][peak_index])
    fwhm_values.append(fwhm)
    #print(f"Peak at index {peak_index}: FWHM = {fwhm}")

for i in range(len(lista_picos_index)):
    largo_pico(fwhm_values[i], df, lista_picos_index[i])

#largo_pico(lista_distancia_entre_picos[1], df,1,1)
#print(lista_picos_index[1])
#print(calculate_fwhm(lista_picos_index[1]))

fig.show(renderer='notebook')

lista_de_angulos_2theta =  angulos_peaks
lista_de_angulos_theta = []
for angulo in lista_de_angulos_2theta:
    lista_de_angulos_theta.append(angulo/2)

sen_angulos_plot = []
cos_angulos_plot = []
#print(beta_constant_value)
#beta_constant_value = 3
contador = 0
for angulo in lista_de_angulos_theta:
    #print(fwhm_values[contador])
    #print(np.deg2rad(angulo))
    # Convert the angle from degrees to radians
    deg_angle_sen = np.sin(np.deg2rad(angulo))
    deg_angle_cos =  np.cos(np.deg2rad(angulo))
    
    sen_angulos_plot.append(4 * deg_angle_sen)
    cos_angulos_plot.append(fwhm_values[contador] * deg_angle_cos)
    #print(np.cos(angulo) * beta_constant_value)
    contador += 1

#print(sen_angulos_plot)

# Add a scatter plot of the list values
# Create a Plotly figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=sen_angulos_plot, y=cos_angulos_plot, mode='markers', marker=dict(color='blue')))

# Update layout if needed (e.g., title, axis labels)
fig.update_layout(title='SrTiO3 ßcosø vs senø', xaxis_title='4senø', yaxis_title='ßcosø')

# Show the plot
fig.show(renderer='notebook')

# Se crean arreglos de numpy con las listas de los angulos senos y cosenos ya tratados
x = np.array(sen_angulos_plot)
y = np.array(cos_angulos_plot)

# Se hace la regresión lineal con numpy
slope, intercept = np.polyfit(x, y, 1)  # 1 for linear regression

# Creamos la función de la linea de regresión con la variable slope e intercept
regression_line = slope * x + intercept

# Creamos una gráfica
fig = go.Figure()

# Graficamos los puntos originales
fig.add_trace(go.Scatter(x=x, y=y, mode='markers', name='Datos originales'))

# Graficamos la linea de regresión
fig.add_trace(go.Scatter(x=x, y=regression_line, mode='lines', name='Linea de regresión'))

# Agregamos las etiquetas
fig.update_layout(
    title='SrTiO3 ßcosø vs senø con Regresion lineal',
    xaxis_title='4senø',
    yaxis_title='ßcosø'
)
equation = f'y = {slope}x + {intercept}'

# Calculamos R cuadrada
r_squared = r2_score(y, regression_line)
# Mostramos la gráfica
fig.show(renderer='notebook')

display(Markdown(f'''
&nbsp;

<div align="center"> Ecuación obtenidad: {f'y = {slope}'}<b> x </b>+ {f'{intercept}'} </b></div>
<div align="center"> $R^2$= {r_squared} </b></div>

&nbsp;

'''))

K = 0.9 # Para cúbicas según la referencia
LAMBDA = 1.5406 # Longitud de onda Cobre K alfa

wh_cristalito = (K * LAMBDA)/intercept

    
display(Markdown(f'''
&nbsp;

<div align="center"> Tamaño de cristalito calculado por Williamson-hall : <b> {wh_cristalito} nm </b></div>

&nbsp;

'''))



if slope > 0:
    display(Markdown(f'La pendinente es positiva por lo tanto podemos argumentar que el sistema tiene tiene una __tensión__'))
else:
    display(Markdown(f'La pendinente es negativa por lo tanto podemos argumentar que el sistema tiene tiene una __compresión__'))

x_values = df[' 2THETA'].values
y_values = df[' PSD'].values


# Se calcula el área bajo la curva (regla trapezoidal)
area_total = np.trapz(y_values, x=x_values)

# Se grafíca la señal del difractograma
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_values, y=y_values, mode='lines', name='XRD original'))

# Creamos un poligono que dibuje el área bajo la curva
x_polygon = np.concatenate([x_values, x_values[::-1]])
y_polygon = np.concatenate([y_values, np.zeros_like(y_values[::-1])])

fig.add_trace(go.Scatter(
    x=x_polygon,
    y=y_polygon,
    fill='tozeroy',
    mode='none',
    fillcolor='rgba(250, 100, 80, 0.5)', 
    name=f'Area bajo la curva: {area_total}'
))

# Customize layout
fig.update_layout(
    title='Difractograma con el área bajo la curva',
    xaxis_title='2 theta',
    yaxis_title='PSD'
)

# Show the plot
fig.show(renderer='notebook')

display(Markdown(f'''
&nbsp;

<div align="center"><h3> Area bajo la curva obtenida: {area_total}</h3></div>


&nbsp;

'''))

# Calculate the area under each peak using np.trapz
peak_areas = []

for i in range(len(peaks) - 1):
    peak_start = peaks[i]  # Se calcula el inicio del pico
    peak_end = peaks[i + 1] if i + 1 < len(peaks) else len(df)  # Se calcula el final del pico
    
    # Se obtienen los valores para x y y
    x_peak = df[' 2THETA'][peak_start:peak_end].values
    y_peak = df[' PSD'][peak_start:peak_end].values
    
    # Calcula el area usando la función np.trapz
    area = np.trapz(y_peak, x=x_peak)
    peak_areas.append(area)

# Imprime el area calculada de cada pico
for idx, area in enumerate(peak_areas, start=1):
    print(f"Area bajo el pico {idx}: {area}")

print(f'\nSuma del area te todos los picos: {sum(peak_areas)}')

Area bajo el pico 1: 9736.148850000009
Area bajo el pico 2: 11060.014950000117
Area bajo el pico 3: 5688.808999999912
Area bajo el pico 4: 4244.044500000017
Area bajo el pico 5: 3367.6622500000008
Area bajo el pico 6: 5740.697700000029
Area bajo el pico 7: 4104.956550000022
Area bajo el pico 8: 1971.0915000000175
Area bajo el pico 9: 1454.2565000000009

Suma del area te todos los picos: 47367.68180000013

porcent_crystallinity = (sum(peak_areas)/area_total) * 100

display(Markdown(f'''
&nbsp;

<div align="center">% de cristalinidad obtenida: <b> {round(porcent_crystallinity, 3)} %</b></div>


&nbsp;

'''))

XRD File reader and plotter¶

Lectura del Archivo Crudo¶

Conversión del Archivo¶

Generación de Gráfica para la Fórmula de Debye-Scherrer¶

Cálculo de tamaño de cristalito a partir de la formula de Debye-Scherrer¶

Generación de Gráfica para la Fórmula de Williamson-Hall¶

Se obtienen los datos necesarios para hacer la gráfica $\beta cos\theta$ vs $sen\theta$¶

Se grafíca $\beta cos\theta$ vs $sen\theta$¶

Regresión lineal de los datos¶

Calculo del cristalito utilizando la fórmula de Williamson-hall¶

Cálculo del Porcentaje de Cristalinidad¶

Area bajo la curva obtenida: 78767.84545

Se obtiene el area bajo la curva pero solo de los picos de la gráfica¶

Cálculo del % cristalinidad¶

Conclusión¶